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Abstract 

We consider a process in which there are p-species of particles, i.e. A\,Ai, • • • , A p , 
on an infinite one-dimensional lattice. Each particle Ai can diffuse to its right neighbor- 
ing site with rate Di, if this site is not already occupied. Also they have the exchange 
interaction Aj + Ai — > Ai + Aj with rate rjj . We study the range of parameters (inter- 
actions) for which the model is integrable. The wavefunctions of this multi-parameter 
family of integrable models are found. We also extend the 2-species model to the case 
in which the particles are able to diffuse to their right or left neighboring sites. 
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1 Introduction 



Our understanding of nonequilibrium statistical physics is far behind that for the equilibrium 
theory. Even simple models may pose a formidable problem if one wants to approach them 
analytically. As an interesting example of stochastic models, which may be investigated 
analytically in some few cases, are one-dimensional reaction-diffusion processes which are of 
both theoretical and experimental interest in a very wide context of physics and chemistry, 
such as stochastic spin flip dynamics traffic flow 0, |3|, the kinetics of bipolymerization 
d, U, reptation of DNA in gels || [7|], interface growth f|, diffusion in zeolites |I| [TTfl , 
and many other phenomena. 

Asymmetric Simple Exclusion Processes (ASEP) in one dimension, is one of the simplest 



example of a driven diffusion system |12|, [13|. For example, the totally ASEP model describes 
a process in which each lattice site can be occupied by at most one particle and the particles 
hop to their right neighboring site if they are not already occupied, with a rate which is 
the same for all particles, otherwise the attempted move is rejected. The dynamics of these 
models can be fully specified by a master equation and an appropriate boundary condition, 
which imposed on the probabilities appear in the master equation. Using the coordinate 
Bethe ansatz, the author of |14| has exactly obtained the iV-particle conditional probabilities 
of totally ASEP, in which the particles can move to left and right with different rates. 

Now the interesting point is that if one changes the boundary condition, without altering 
the master equation, one can model another reaction-diffusion processes even with long 



range interactions. For example in [|TjJ, the so called "generalized totally ASEP model" 
has been exactly solved in this way. In this model the particle hops to the next right site 
by pushing all the neighboring particles to their next right sites, with a rate depending on 
the number of right neighboring particles. The partially generalized ASEP model has also 



been studied in [16]. Note that in all of these cases, the solvability of the models is shown 
by proving the factorization of TV-particle S-matrices into 2-particle ones, which the latter 
were found exactly. 

In all the above ASEPs, there is only one species of particle, that is all the particles are of 
the same type. But if one considers the two, or more, species problems, the situation becomes 
more complicated. The main complexity arises from the fact that the above mentioned 
factorization of iV-particle systems, reduces to the condition of satisfying the two-particle 
S-matrices in the Quantum Yang-Baxter Equation (QYBE). In one species models, the 
S-matrices are not really matrices, they are c-numbers and therefore their satisfying in 
QYBE becomes trivial. This new condition can hardly restrict the number of solvable 
models with more than one-species particles. In [17[] , a class of two species reaction-diffusion 
processes with following properties has been considered: 1) the particles diffuse to their right 
neighboring sites, 2) they can be annihilated or created, but the total number of particles 
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are constant, and 3) the interaction rates are all the same. It is shown that among 4096 
types of models with the above properties, which can be modeled by a master equation and 
a number of boundary conditions, there are only 28 independent interactions which their 
two-particle S-matrices satisfy the QYBE and therefore are solvable. The third condition 
(equality of the interaction rates) was very crucial in proof of solvability. 

In this paper we want to study the effect of interaction rates in solvability of p-species 
reaction-diffusion processes, by considering a specific model. We begin our investigation 
by choosing one of the two-species interactions that has been introduced in ]T7|, but with 
different interaction rates, and try to obtain the range of parameters to insure the solvability 
of the corresponding extended p-species model. As we show, we must restrict ourselves to 
a narrower and narrower range of parameters, as we go ahead, and finally arrive at a model 
with a specific relation between the interaction rates and also a specific range for these rates. 

The plan of the paper is as following. In section 2, we begin with following two-species 
reaction-diffusion processes: 

A + 0^40 + A, 
5 + 0^0 + B, 

A + B^B + A, (1) 
B + A^ A + B, 



and write down a master equation and a number of boundary conditions to describe the 
dynamics of these interactions. Da and Db are the right-diffusion rates of A and B particles, 
respectively, and s and r are the rates of transforming (exchanging) A and B particles to 
each other for (...AB...) and (...BA...) configurations, respectively. We will show that only 
for Da = Db = 1 case, there exists the coordinate Bethe-ansatz solution for probabilities 
(note that taking D = 1 is in fact a choosing of time scale). Moreover, we will see that the 
consistency of the solutions (which will be appeared as satisfaction of 2-particle S-matrix in 
QYBE) restricts us to r = or s = cases (which are the same after relabeling A < — > B). 
We therefore conclude that the solvable model (interaction) is the following one-parameter 
family process 

A + <ft±$ + A, 

5 + 0^0 + 5, (2) 
B + A^ A + B. 

Note that for r = 1, interactions (0) are one of the 28 interactions introduced in [17]. We then 
generalize the interactions (fj) to a p-species model in which the particles Ai(i = 1, ■■■,p) can 
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diffuse to their right neighboring sites, all with equal rate one, and also they have exchange 
interactions with different rates: 

A + ^ + A, 

Aj+Ai^Ai + Aj, (j>i). (3) 

We label the species such that the configuration (...A, A...) can g° to (.. AjA^-...), only when 
j > i, therefore r in eq.(fj) is in fact r\ 2 . After a lengthy calculation, we show that there 
must be a specific relations between r^s until the p 2 x p 2 two-particle S-matrix satisfying 
the QYBE. 

In section 3, we calculate the two-particle conditional probabilities of reaction (0), and 
show that only for < r < 2 we are able to calculate these probabilities by a standard 
superposition of the eigenfunctions with real eigenvalues. The long-time behavior of the 
probabilities is also discussed. Finally in section 4, we generalize the reactions (0) to the 
case where the particles can diffuse to both right and left, as following: 

A + 0^0 + A 
B + Q^Q + B, 

+A 2 A+ 0, (4) 
0+5^5+ 0, 
B + A-^ A + B, 
A + B A B + A. 

We show that there must be a fine tuning of parameters if one demands the reactions (4) to 
be solvable. 



2 p— species exchange— diffusion processes 
2.1 The master equation for 2— species case 

Consider the interactions introduced in eq.(|l]). The basic quantities that must be calculated 
are the probabilities P ai a 2 ...a N (xi, x 2 , x^] t) for finding at time t the particle of type a.\ at 
site x\, particle of type a 2 at site x 2 , etc.. Each a.i can be A or B. Following ||14|| , we take 



these functions to define probabilities only in the physical region x\ < x 2 < ... < xn, and 
the regions where any two adjacent coordinates are equal, are the boundaries of the physical 
region. For — Xi > l,Vi, the particles can only hop to their right neighboring sites and 
therefore the master equation is: 

d N 
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t)l (5) 

where the first N terms are the sources of P ai a 2 ...a N ( x i, x 2, ■■, %n] t) and the second N terms 
are the sinks of it. It is obvious that if Xj+i = Xi + 1 for some z's, then some of the probability 
functions in the right hand side of eq.(|J) go out from the physical region. So we need to 
specify the boundary terms. The specification of these terms depends on the details of 
the interactions of the particles. For exchange interactions defined in eq. ([!]), the suitable 
boundary conditions are: 

D A P BA (x, x) = sP AB (x,x + 1) + (D B - r)P BA (x,x + 1), 
D B P AB (x,x) = rP BA (x,x + 1) + (D A - s)P AB (x,x + 1), (6) 
P aa (x,x) = P aa {x, x + 1) , (a — A,B), 

in which the time variable and all the other coordinates have been suppressed for simplicity. 
To justify these boundary conditions, it is enough to examine them in some specific cases. 
Let us do it for a rather complicated case, for example P ABBA (x, x + l,x + 2, x + 3). From 
master equation ([5]), we have, 

d 

tt;Pabba{x, x + 1, x + 2, x + 3) = D A P ABBA (x - 1, x + 1, x + 2, x + 3) 
ot 

+ D B P ABBA (x, x, x + 2, x + 3) + D B P ABBA (x, x+l,x + l,x + 3) (7) 

+D A P ABBA (x, x + l,x + 2,x + 2)- 2(D A + D B )P ABBA (x, x + 1, x + 2, x + 3). 

If we use the relations (6) in the second, third and fourth terms of the right hand side of (7), 
we find 

d 

-^Pabba(x, x + 1, x + 2, x + 3) = D A P ABBA (x - 1, x + 1, x + 2, x + 3) + 

rP BABA (x, x+1, x+2, x+3)+sP ABAB (x, x+1, x+2, x+3)-(D A +r+s)P ABBA (x, x+1, x+2, x+3). 

(8) 

This equation is exactly what we expect from interactions (|I|), because the source terms of 
configuration (.J ABB A®...) are: (...AflBBA®...) (with rate D A ), (.JBABA®...) (with 
rate r) and (...®ABAB$...) (with rate s), and its sink terms are: (...$ABB$A...) (with rate 
D A ), (...0SASA0...) (with rate s) and (...0ABAB0...) (with rate r). It can be shown that 
the boundary conditions (6) results the correct terms for any desired configuration. 

2.2 the Bethe ansatz solution (2— species) 

Now we want to solve the master equation (5) with boundary conditions (6) by the coordinate 
Bethe ansatz method. First we define ^ ai ...a N (xi, ---^tv) through, 

P ai ...a N {xu %N] t) = e" eiV '^ ai ... aiV (xi, x N ), (9) 
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and then substitute it in eq.(5), which results 



JV 



N 



J2 D a l ^a 1 ...a N (Xl, Xi-i, Zj - l,X i+1 , ...,X N ) = D ai - € N )^ ai ... aif (xi , ...,X N ). (10) 



i=l 



To solve this equation, we use the coordinate Bethe ansatz for each of the components 



(xi, ...,Xat) : 



'a) = £4"— V^)- x , (11) 



where x and p stands for A-tuples coordinates and momenta, respectively, and cr(p) is a 
permutation of momenta. The sum is over all permutations. Inserting (11) into (10) yields: 



E 



N N 
3 = 1 .7 = 1 



^(ai...ajv) e icr(p).x _ q 



'121 



As ^K-«iv) e Mp)-x s 

are linearly independent for different a's, the only solution of eq.(12) is: 

AT N 



e N -J2D aj +E D aj e~^ = 0, V<r, 
3=1 j=i 



(13) 



or 



N 



N 



N 



D aj e- 1 ^ = D aj e~ 1 ^ = - = E A^e"""^, (14) 
j=i j=i j=i 

where n is the number of elements of permutation group. As these equalities must hold for 
an arbitrary p, the only nontrivial solution is: 



D A = D B = 1. 



(15) 



Now as for any group element a we have 

TV 



J- e -™M = e ~ ipi + ... + e~ ipN , 

3=1 



(16) 



the equalities (14) satisfy satisfactorily. Therefore, the Bethe ansatz solution exists only for 
equal diffusion rates, which in this case the eigenvalue is found to be (by eq.(13)): 



N 



€ n — E(i — e~ ipj )- 

3=1 



(17) 



It is easier to consider i & ai ... aN (xi, —iXn) and A^* 1 ~ aN ' > as the components of the tensors * 
and Ao-with rank N, respectively, therefore eq.(ll) can be written as 



*(£!,..., x N ) = ^A CT e^ p )- x . 



(18) 
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The boundary conditions of \I> can be obtained by substituting eq.(9) in (6), with D A = 
D B = 1. The resulting equation is 



C,C,-) = Kk+i*(->CC + i,-) 5 



(19) 



where 



with 



>k,k+l 



k,k + l 



(20) 



(21) 



/ 1 \ 

1-s r 

s 1 -r 

V o o o i y 

The coefficients A CT 's in eq.(18) must be found by substituting the wavefunction (18) into 
the boundary condition (19), which yields 

J2 eT (Pj) x j+ i ( eT (Pk)+v(Pk+l))t 



'1 — e ia ( ph+1 ^h 



0. 



(22) 



As the exponential part of eq.(22) is symmetric with respect to p k < — > Pk+i, if we a ls° 
symmetrize the remaining terms with respect to this interchange, we obtain 



(1 _ e ^^b kM1 )A a + (1 - e^b kM1 )A aak = 0, 



(23) 



where a k represent the permutation group element which only interchanges p k and Pk+i- 
Therefore 



where 



in which 



S ktk+1 {a{p k ), a(p k+1 ))A a , 

h,k+i(zi,z 2 ) = 1 <S> ■ • • <8> 1<S> S(z^, z 2 ) ®1( 

k,k + l 

S(z 1 ,z 2 ) = -(l-z 1 b)-\l-z 2 b). 



(24) 
(25) 
(26) 



In above equations, z k stands for e tpk . In this way, all the A CT coefficients are determined 
in terms of Ai which is fixed by the initial conditions (the particles' positions at t — 0). It 
seems that we have solved the problem for arbitrary b (i.e. interaction), but it is not true 
(note that we have not yet used the explicit form of b in deriving A a ). The crucial point 
is that o\o 2 (J\ and a 2 a±a 2 are equal as elements of permutation group, therefore we should 
impose the following condition on the corresponding A CT 's: 



A = A 

Jr% -U\(72(J\ - r *-0"20"l(T2 7 



(27) 
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and this highly restricts the allowed b matrices (i.e. interactions ). It can be shown that 
eq.(27) reduces to following relation for S{z ll z 2 ) matrices (see [|l7j for more details), 

(S(w, t) <g> l){l®S{z, t))(S(z, w) ® 1) = (1 <g> S(z, w))(S(z, t) <g> 1)(1 ® S{w, t)), (28) 

in which z = e tpi ,w = e tp2 , and t = e m . Note that eq.(28) is nothing but the Quantum 
Yang-Baxter equation. 

Now if one calculates the S-matrix from eq.(26), using b from (21), the QYBE (28) 
reduces to a 8 x 8 matrix with fourteen nonzero elements (after writing eq.(28) as RHS - 
LHS = 0) that must be equated to the zero matrix. The elements are functions of pi,P2,P3, f 
and s which all must be equal to zero for arbitrary momentum values Pi,P2 and p%. It can 
be shown that the unique solutions of these fourteen equations are: 

solution 1 : r = 0, arbitrary s. 

solution 2 : s — 0, arbitrary r. (29) 

As these two solutions are equivalent (by relabeling A < — > B), so the only integrable model 
is the one indicated in eq.(2). 

2.3 The p— species model 

Now let us generalize the reaction (2) to the case where there exists p kinds of particles, which 
we label them by A±, A 2 , A p . Each particle can diffuse to its right neighboring site, and any 
two particle can exchange with each other. From the results obtained in previous subsection, 
we know that if we want the model to be integrable, we must restrict ourselves to the case 
where the particles' diffusion rates are equal (scaled to one), and also for each two particles 
there is only one allowed exchange interaction. For example A 2 + A\ — > Ax + A 2 is allowed, 
but Ai + A 2 — > A 2 + Ai is forbidden. We label the particles such that Aj + A{ — > A{ + Aj is 
allowed only for j > i, and denote the reaction rate of this interaction by r^, see eq.(3). 

The master equation of this p-species model is again eq.(5) (with D a . = 1), but now 
each ai can be A±, A 2 , and A p . The boundary conditions are the generalization of eq.(6) 
but now for each two particle species Ai and Aj, i.e. Pab — * Pij, Pba — * Pji, s = 0, 
Da = P>b = 1, and r — >■ r^-; so 

Pij{x, x) = Pij{x, X +1) + TijPji{x, X + 1), j > i, 

Pji(x,x) = (l-rij)Pji(x,x + l), j > i, (30) 
Pa(jx, X^j Pa^X , X ~t~ 1). 

The wavefunctions can be again factorized by the eq.(9), and the coordinate Bethe ansatz 
solution (18) is yet valid. The boundary conditions can be rewritten as (19), but here b is 
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the following p 2 x p 2 matrix: 

b = E u <g> Ejj + J2 r, J F., l ® /•;,-,• + £(1 - r^E u <g> (31) 

i<j i<j i>jr' 

where EV,- is a p x jo matrix with elements (Eij)ki = Sik^ji- It can be shown that the S-matrix 
(26) becomes: 

v l-iyfl-r'-) p z-w 

= ,£ (1-4)^-1 ^ 8 s " + .£ 8 E * (32) 

in which, 

r H={ ■( ■ / ■ • ( 33 ) 
[ nj if i < j 

Now expression (32) must satisfy the QYBE (28). After a lengthy calculation, it can be 
shown that the only nontrivial solutions of QYBE are as following (for each i < j < k 
indices): 



solution 1 
solution 2 
solution 3 



Tij = 0, Tik and Tjk arbitrary, 

Tik = 0, rjk = 0, and arbitrary, (34) 
Tij = r ik , r jk arbitrary. 



For any set of interaction rates where each three of them satisfy any of the solutions 1, 2 or 
3, with the constraint that the relations between r^s must be consistent in all subsets, we 
have an integrable p-species model with wavefunction (18) whose coefficients are determined 
by eq.(24) and S-matrix introduced in (32). 

For p — 3, the allowed sets of interaction rates are exactly the same as the three solutions 
(34) with (ijk) = (123). But for p > 3 cases, we can choose different consistent solutions 
for any (ijkys and therefore extracting all the allowed sets are not so easy. For example for 
p = 4, in which there are six interaction rates r 12 , r 13 , r 14 , r 23 , r 24 , and r 34 , the allowed sets 
of parameters are as following: 

{^14,^24,^34} , 

{ri 2 ,r 34 } , 
0~i 3 , r 23 } , 

{^13,^24} , 

{ri4,r 23 } , 
{ri4,r 3 4,r 2 3 = r 24 } , 
j> 24 ,r 34 ,ri 2 = r i4 } , 

{r 2 4,r 34 ,ri 3 = r u } , (35) 
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{^24,^34,^12 = r 13 = r lA ] , 

{^13,^23 = ^24} , 

{r 23 , r 12 = r 13 } , 

{^23,^13 = ^14} , 

{r 2 3,n 2 = = r u } , 

{^34,^13 = ^14,^23 = r 24 } , 

{r 3 4, ri2 = ri3 = r M , r 23 = r 24 } . 

Note that in all the above allowed sets, we have only brought the free parameters and the 
relations that must be satisfied by them, and the zero reaction rates have not been written. 
In this way we find a large class of multi-parameter p-species integrable reaction-diffusion 
models. 



3 two— particle conditional probabilities for 2— species 
model 

Now for the simplest case, that is the 2-species reactions (2), let us calculate the two-particle 
conditional probabilities P(a\, a 2 , Xi,x 2 ; t\P%, /3 2 , 2/1,1/2; 0), which is the probability of finding 
particles a.\ and a 2 at time t at sites X\ and x 2 , respectively, if at t = we have the particles 
/3i and f3 2 & t sites yi and y 2 , respectively. These probabilities can be found by a linear 
combination of eigenf unctions P(xi,x 2 ). Therefore, 



(x;t|/3,y;0) 



/ Paa \ 
Pab 
Pba 

V Pbb ) 

I /(pi,p 2 )e" e2 '*(xi, x 2 )dp 1 dp 2 



±v f e ~ e2t e- ip y 



(36) 



(27T) 





( a ) 




( a \ 


< 


b 


e i( Pl X 1 +P2X 2 ) _|_ S 12 {pi,p 2 ) 


b 




c 




c 








[dj 



> dpidp 2 . 



In these expansion, P(x; t\@, y; 0) stands for P(«i, a 2 , Xi, x 2 \ (3 2 , yi, y 2) 0) and f(pi,p 2 ) is 
the coefficient of expansion, where in the second equality we choose it to be -^p- / e~ ip y (see 



14j - |17|1 ). e 2 = 2 — e ipi — e ip2 (see (17)) and is the two-particle wave function (18), in 
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which eq.(24) has been used for A^A^ = £12 (pi, p2)Ai). The column matrix 



/ a \ 
b 

c 



stands 



for Ai, which its components must be determined by initial condition and Si2(pi,P2) is: 

/ si \ 



Sl2(Pl,P2) 



si s 2 
s 3 
\ si J 



(37) 



where 



«1 



S2 



S3 



e ipi _ 1 

r(e ip2 - e ipi ) 
(1 - e^ 1 ) [1 + (r - l)e*Pi]' 
(1 - r)e ip2 - 1 
1 + (r - l)e*Pi ' 



(38) 



The matrix <Si2(pi,P2) is obtained from eq.(26) in which the matrix b in eq.(21) (with s = 0) 
has been used. By inserting eq.(37) into eq.(36), we find: 



/ Paa \ 
Pab 
Pba 
V Pbb ) 



(x;t|/?,y,0) 



/ a(F (t) + F 1 (t)) 
b(F (t) + F^t)) + cF 2 (t) 
c(F (t) + F 3 (t)) 
d(F (t) + F 1 (t)) 



(39) 



in which, 



Fo(t) = 
Fi(t) = 



(2tt)= 
1 



J e -^ e ip-(x-y)d Pl d P2 , 



e -^ e i$*-py) Si ( Pl , P2 ) dpidp2 , (i = i )2 ,3). 



(40) 
(41) 



In above equations, we have suppressed the x and y dependence of Fj's, for simplicity, and 
P= (p 2 ,pi). Now at t — 0, the configuration of the system can be one of the (A, A), (A, B), 
(B, A) or (B, B), where the first particle is at site y\ and the second one at y 2 , therefore the 
only acceptable behavior of i^(0) (i = 0, 1, 2, 3) are: 



F o (0) 

Fi(0) = F 2 (0) = F 3 (0) = 0. 



(42) 



Fo(0) is obviously correct (see eq.(40)). For other Fj's, first we must set P \ — > pi + to 
avoid the singularity arising from e ipi — 1 term in denominator of S\ and S2 (see [14-17]). In 
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this way one can show that Pi(0) = 0. But F 2 and P 3 have another singularity because of 
1 + (r — l)e tpi term in denominator of s 2 and s 3 . One can easily show that this singularity 
can be avoided only when 

< r < 2, (43) 

and for this range of interaction rates, we have F 2 (0) = P 3 (0) = 0. Therefore the validity of 
expansion (36) is restricted to range (43). At t ^ 0, we find 



Fo(t) 
Fi(t) 



-2t 



^xi-yi 



^xi—yi 



-It 



xi - yi)\ (x 2 - y 2 )!' 



x 1 -y 2 + l)\ (x 1 -y 2 )\ 
1 



oo -i-X2-yi+k 

y - 



fxx-y 2 oo 

W) = re- 2 y £ 

{xi-y 2 )\ l j^ Q 



x\ — y 2 + 1 x 2 - yi + k + / + 1 



(l-r) 1 



(44) 

£X2-yi+k+l+l 

(x 2 - yi + k + l)V 



F 3 (t) = e 



-2t 



£X!-y2 

(xi - y 2 )! 



r)t 



(x 1 -y 2 + l)\ 



00 +X2~yi+k 

Y(l-r) k - 

to (x 2 -y 1 + k)\ 



One can now obtain the two-particle conditional probabilities for different initial condi- 
tions: 

1. If at t — 0, the particles j3i = (3 2 = A were at y\ and y 2 , respectively, we must take 
a = 1 and b = c = d = 0. So at t 7^ 0, we have 



P A A(x;t\A,A,y;0) = F o (t) + F 1 (t), 



(45) 



and all other P's are zero. 

2. If Pi = A and f3 2 = B, we must take b = 1 and a = c = <i = 0. Therefore the only 
nonzero probability is 

P AB (x; t\A, B, y; 0) = F (t) + F 1 (t). (46) 

3. If Pi — B and fi 2 = A, we have c = 1 and a = b = d = 0. So = Pb B = and 



P AB (x;t\B,A,y;0) 
P BA (x;t\B,A,y;0) 



F (t)+F 3 (t). 



(47) 



4. And finally if (3i = (3 2 = B, we have d — 1 and a = 6 = c = 0. So the only nonzero 
probability is 

P BB (x; i|P, P, y; 0) = F Q (t) + F^t). (48) 

Note that the appearance of the above probabilities is in agreement with our reactions (2). 

As another check of our results, it may be interesting to study the long time behavior of 
these probabilities, in special P4#(x; t\B, A, y; 0) which is the only nondiagonal nontrivial 
case. We expect that if at t = we have a B particle at site yi and an A particle at site y 2 
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(with 2/2 > Hi), we must certainly have two B particles at t — > oo somewhere at y± < X\ < x 2 
and y 2 < £ 2 < 00 sites. In other words, we expect 

f; £ P j4B (x;^oo| J B,Ay;0)^l. (49) 

After some calculations, one can show that 

OO X2 — 1 OO 

E E P ab(x; t|fl, A, y; 0) = e" 2 ' ]T [l - (1 - r)^ 1 ] (2t) + 7 n+2/2 _ 2/1+1 (2t)] , 

i2=!/2 ii=!/i n=0 

(50) 

where I„(x) is the n-th order Bessel function of the first kind. To obtain the long-time 
behavior of (50), one may simply use the following asymptotic form of I n (x) at x — > 00, 

Uz) "> -j=, (51) 

in eq.(50). But it is not correct since eq.(51) is only valid for x > n, but in eq.(50) we have 
a sum over n where for every large definite t, there exist infinite number of n which are 
greater than t (it can be shown that if one calculates this limit without noting this point, 

finds infinity for eq.(50), which is obviously wrong). To calculate this limit, if one uses the 

00 

identity E I n {x) = e x and takes advantage of equality I n (x) = I- n (x), can show that 

n=— 00 

00 k— 1 

E I n +k(x) = \[e x - E I n {x)). Therefore 

n=0 " n=—k+l 



e 



00 1 k— 1 00 

-» £ [l - (1 - r)" +1 ] I n+k (2t) = --e~ 2t J2 Wt)-e- 2t J2(l-r) n+1 I n+ k(2t). (52) 
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ra=0 ^ ra=-fc+l n=0 



Now in the second term of the RHS of (52) n is bounded, so eq.(51) can be used for it which 
leads to zero in t — > 00 limit. For the third term we note that — 1 < 1 — r < 1 (see eq.(43)), 
so (1 — r) n+1 — > for large n . One can show that this extra (1 — r) n+1 factor causes the 
third term in RHS of (52) goes also to zero at t — > 00 limit. Therefore the t — > 00 limit of 
RHS of (51) is equal to 1/2, from which eq.(49) is proved. 

4 2-species model with left-right diffusion 

In this section we want to study the range of parameters (reaction rates), which makes the 
reactions (4) integrable. In this case the master equation is 

d N 

at i=l 

JV 

+-D_L ^ ^"ai...ajv( a; l; •••) ^i-lj + 1) x i+l, •••) ^TVj ^) 
i=l 

t), (53) 
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in which we have used a time scale so that 



D R + D L = 1. (54) 

The boundary conditions are 

D R P AB {x,x) + D L P AB {x + l,x + 1) = rP BA {x,x + 1) + {1 - s)P AB {x,x + 1), 
D R P BA (x,x) + D L P BA (x + l,x + l) = sP AB {x,x + l) + {l-r)P BA {x,x + l), (55) 
D R P aa (x,x) + D L P aa (x + l,x + I) = P aa {x,x + 1), (a = A,B). 

Note that for D L = 0, eqs.(53) and (55) reduce to eqs.(5) and (6), respectively, and for 
r = s = 0, these equations lead to corresponding ones in []14| . 

We must proceed the same steps as previous section, where in this case lead to the 
following relations for energy and boundary conditions: 

N 

e N = Y,(l-D R e- i *-D L e i n), (56) 
i=i 

and 

D R V{...,(, C, -..) + D L *(...,C + 1, C + 1, •••) = b M+1 *(...,C, C + 1, •••), (57) 

where bfc^+i is defined through (20) with b in eq.(21). If we substitute the coordinate Bethe 
ansatz (18) in eq.(57), the relation between coefficients is like eq.(24), but now with following 
S-matrix: 

S(zi, z 2 ) = -{D R + z 1 z 2 D L -z 1 h)-\D R + z lZ2 D L -z 2 b). (58) 

This S-matrix must satisfy the QYBE (28). Like the previous case, here there are also 
fourteen equations that must be solved for r, s and an extra D R parameters {Dl is fixed 
by D R through eq.(54)). These equations are highly nonlinear and we are not able to solve 
them exactly, even by using the standard programs like MAPLE. So we restrict ourselves to 
the cases in which 

r + s = l. (59) 

In this way, we can find the complete set of solutions. We believe that there are no other 
solutions even if the constraint (59) is removed (we have checked many other cases, but no 
one satisfied QYBE). The solutions with r = or s = are not new. They are the known 



models like the models introduced in (29), or the simple diffusion models introduced in [H] 
(in which both r and s are zero). There are only two new integrable models, as follows: 

B + ^ + B, 

0, (60) 
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d}+B^B+ 0, 
B + A^A + B, 
A + B ^ B + A, 

and 

£ + 0^0 + 5, 

(61) 

®+B^B+ 0, 
B + A^A + B, 
A + B^ B + A, 

in which Dr = 1 — Dl- In this way we find two one-parameter family integrable models, 
which their time dependent probabilities can be found by eq.(9), with in eq.(56), ^ is 
given by (18) and (24), with S-matrix introduced in (58), and b in (21). 
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